#include "dyn_multifunctor_multiMaps_ikedaMap.h"

IkedaMap::IkedaMap(const double R, const double C_1, const double C_2, const double C_3, const int p)
    : MultiMaps(p)
{
    setR(R);
    setC_1(C_1);
    setC_2(C_2);
    setC_3(C_3);
}

colvec IkedaMap::iterate(const colvec X,int p)
{
    if (p==0)
        return X;

    // We declare some variables that will be useful
    // in the computation.
    double tau;
    double x_temp;
    double y_temp;

    // This will be the return vector.
    colvec answer(2);
    answer = X;

    for (int i=0;i<p;i++)
    {
        tau = c_1 - c_3/(1+accu(pow(answer,2)));
        x_temp = r+c_2*(answer(0)*cos(tau)-answer(1)*sin(tau));
        y_temp = c_2*(answer(0)*sin(tau)+answer(1)*cos(tau));

        answer(0) = x_temp;
        answer(1) = y_temp;
    }
    return answer;
}

mat IkedaMap::jacobian(const colvec X, int p)
{
    /*! The jacobian of a pth iterate is given
     * by the multiplication of the jacobians
     * of the p iterates of the point.
     * We therefore generate a matrix of the p
     * iterates of the point. We then
     * calculate the jacobian at these points
     * and multiply the resulting matrices.
     */

    // Used to hold the p iterates temporarily.
    vec temp;

    // Return matrix and parts of it.
    mat jac(2,2),t_jac(2,2), part1(2,2), part2(2,2);

    for (int i=0;i<p;i++)
    {
        // We calculate the ith iterate.
        temp = iterate(X,i);

        // We fill in the jacobian matrix for this point.
        double tau = c_1 - c_3/(1+accu(pow(temp,2)));

        // We break it up in two for readability.
        part1 << 1 << -tan(tau) << endr
              << tan(tau) << 1 << endr;

        part2 << temp(0)*(temp(0)*tan(tau)+temp(1)) << temp(1)*(temp(0)*tan(tau)+temp(1)) << endr
              << -temp(0)*(temp(0)-temp(1)*tan(tau))<< -temp(1)*(temp(0)-temp(1)*tan(tau))<< endr;

        t_jac = c_2*cos(tau)*(part1-(2*c_3/(pow(1+accu(pow(temp,2)),2))*part2));

        if (i==0)
            jac = t_jac;
        else
            jac *= t_jac;

    }
    return jac;
}
